# Replication package for 
# "The Economic Leverage of International Organizations in Interstate Disputes"
# Johannes Karreth
# June 30, 2017
# jkarreth@ursinus.edu

# This file: 12_claims_figure.R
# Purpose: Reproduce Figure 1 in the article

# Note: must have read in data before, using 11_claims_readdata.R

# setwd("...")

library("ggplot2")
library("rstanarm")
library("xtable")
library("reshape2")

# Source in functions
source("Functions/MCMClogit.fd.mat.R")
source("Functions/theme_jk.R")
source("Functions/plotBMA.R")
source("Functions/MCMC_roc_prc.r")
source("Functions/geom_flat_violin.R")

# Unit of analysis: claims

# Model 1 (for Figure 1)

m1_eq <- useforce ~ igo_lev3_count_use_bc + salint + saltan + terriss + jointpol7 + rivalry_th + absidealdiff + atopally
m1_mf <- model.frame(m1_eq, data = dat)
m1 <- glm(m1_eq, data = m1_mf, family = binomial(link = "logit"))

m1_stanglm <- stan_glm(formula = m1_eq,
               data = (m1_mf),
               family = binomial(link = "logit"),
               prior = normal(location = 0, scale = 10), 
               prior_intercept = normal(0, 10),
               seed = 20170207,
               cores = 4)

saveRDS(m1_stanglm, file = "Output_MCMC/claims_m1_stanglm.RDS")

# Table

m1_tab <- m1_stanglm$stan_summary[, c(1, 3)]
m1_tab <- data.frame(m1_tab[c(2:9, 1, 11), ])
m1_tab$varname <- c("IGOs with high leverage", 
    "Intangible salience of claim", 
    "Tangible salience of claim", 
    "Territorial claim", 
    "Joint democracy", 
    "Strategic rivalry",
    "UNGA ideal point difference",
    "Allies",
    "Intercept",
    "Log-posterior density")
names(m1_tab) <- c("Mean", "Standard deviation", "Variable")
m1_tab <- m1_tab[, c(3, 1, 2)]

m1_xtab <- xtable(m1_tab, digits = 2)
print(m1_xtab, file = "Output_Tables-and-Figures/claims_m1_table.tex", include.rownames = FALSE)

# First differences (Figure 1)

m1_mm <- model.matrix(m1_eq, data = dat)

temp_sims <- as.matrix(m1_stanglm)
m1.fd <- MCMClogit.fd.mat(model_matrix = m1_mm, mcmc_out = temp_sims, 
                           credint = c(0.025, 0.975), 
                           percentiles = c(0.1, 0.9),
                           full_sims = TRUE)

m1.changes <- matrix(NA, ncol = 2, nrow = ncol(m1_mm))
for (i in 2:ncol(m1_mm)){
  m1.changes[i, ] <- ifelse(length(unique(m1_mm[, i])) == 2 & range(m1_mm[, i]) == c(0, 1), c(0, 1), 
                            quantile(m1_mm[, i], probs = c(0.1, 0.9)))
}

# > m1.changes
# [,1]   [,2]
# [1,]     NA     NA
# [2,] 0.0000 4.0000
# [3,] 1.0000 2.3000
# [4,] 1.0000 5.0000
# [5,] 0.0000 1.0000
# [6,] 0.0000 1.0000
# [7,] 0.0000 1.0000
# [8,] 0.1496 2.6827
# [9,] 0.0000 1.0000

fd.dat <- data.frame(m1.fd)
fd.dat$variable_label <- c("IGOs with high leverage\n(from 0 to 4)", 
                           "Intangible salience of claim\n(from 1 to 2.3)", 
                           "Tangible salience of claim\n(from 1 to 5)", 
                           "Territorial claim\n(no to yes)", 
                           "Joint democracy\n(no to yes)", 
                           "Strategic rivalry\n(no to yes)",
                           "UNGA ideal point difference\n(0.1 to 2.7)",
                           "Allies \n(no to yes)")


fd_long <- melt(m1.fd)

fd_long$Varpos <- ifelse(fd_long$Var2 == "igo_lev3_count_use_bc", 1,
                         ifelse(fd_long$Var2 == "salint", 2, 
                                ifelse(fd_long$Var2 == "saltan", 3,
                                       ifelse(fd_long$Var2 == "terriss", 4,
                                              ifelse(fd_long$Var2 == "jointpol7", 5,
                                                     ifelse(fd_long$Var2 == "rivalry_th", 6, 
                                                            ifelse(fd_long$Var2 == "absidealdiff", 7,
                                                                   ifelse(fd_long$Var2 == "atopally", 8, NA))))))))

fd_long$variable_label <- ifelse(fd_long$Var2 == "igo_lev3_count_use_bc", "IGOs with high leverage\n(from 0 to 4)",
                                 ifelse(fd_long$Var2 == "salint", "Intangible salience of claim\n(from 1 to 2.3)", 
                                        ifelse(fd_long$Var2 == "saltan", "Tangible salience of claim\n(from 1 to 5)",
                                               ifelse(fd_long$Var2 == "terriss", "Territorial claim\n(no to yes)",
                                                      ifelse(fd_long$Var2 == "jointpol7", "Joint democracy\n(no to yes)",
                                                             ifelse(fd_long$Var2 == "rivalry_th", "Strategic rivalry\n(no to yes)", 
                                                                    ifelse(fd_long$Var2 == "absidealdiff", "UNGA ideal point difference\n(0.1 to 2.7)",
                                                                           ifelse(fd_long$Var2 == "atopally", "Allies \n(no to yes)", NA))))))))

# ROPE: use SE of the ratio of 1s to 0s

r_n <- length(m1_mf$useforce)
r_p <- sum(m1_mf$useforce) / r_n
r_se <- sqrt((r_p * (1 - r_p)) / r_n)
ROPE <- r_se

fd.p <- ggplot(data = fd_long, aes(x = reorder(variable_label, -Varpos), y = value)) + 
  scale_y_continuous(labels = function(x) x*100) +
  geom_flat_violin(color = NA, fill = "black") +
  coord_flip() + 
  theme_classic() + 
  ylab("Percentage point change in the probability of using force\nduring a claim as each predictor changes as indicated") + 
  xlab("") + 
  theme(axis.ticks = element_blank(), axis.text = element_text(color = "black"), axis.line = element_blank()) + 
  annotate(geom = "text", y = -0.05, x = 0.15, label = "Use of force less likely", hjust = 1, size = 3) +
  annotate(geom = "text", y = 0.05, x = 0.15, label = "Use of force more likely", hjust = 0, size = 3) +
  annotate(geom = "segment", y = 0, yend = -0.7, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "segment", y = 0, yend = 0.8, x = 0, xend = 0, arrow = arrow(ends = "last", type = "closed"), size = 1) +
  annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray", color = NA, alpha = 0.5) +
  geom_segment(y = 0, yend = 0, x = 0, xend = 10.5, color = "black", linetype = "solid", size = 0.33) + 
  annotate(geom = "text", y = 0.1, x = 7.25, hjust = 0, vjust = 0.5, label = "Region of practical\nequivalence (ROPE)\n[-3, 3]", color = "darkgray", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = ROPE, x = 7.25, xend = 7.3, size = 0.25, color = "darkgray") +
  annotate(geom = "text", y = 0.1, x = 8.25, hjust = 0, vjust = 0.5, label = "% of posterior distribution\nto the left of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.1, yend = -0.2, x = 8.25, xend = 8.15, size = 0.25) +
  annotate(geom = "text", y = 0.35, x = 3.5, hjust = 0.5, vjust = 0, label = "% of posterior distribution\nto the right of ROPE", size = 3) +
  annotate(geom = "segment", y = 0.3, yend = 0.2, x = 3.45, xend = 3.15, size = 0.25) +
  annotate(geom = "text", y = -0.5, x = 7.5, hjust = 0.5, vjust = 0, label = "Mean estimated change", size = 3) +
  annotate(geom = "segment", y = -0.5, yend = -0.275, x = 7.65, xend = 7.9, size = 0.25)

# % < 0

m1.fd_out0 <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0) / length(x), sum(x > 0) / length(x)))
m1.fd_outROPE <- apply(m1.fd, 2, function(x) ifelse(median(x) < 0, sum(x < 0 - ROPE) / length(x), sum(x > 0 + ROPE) / length(x)))

fd_changelab <- ifelse(apply(m1.fd, 2, median) < 0, "Decrease by", "Increase by")
fd_changelab_short <- ifelse(apply(m1.fd, 2, median) < 0, "-", "+")
fd_changetext <- paste(fd_changelab, abs(round(apply(m1.fd, 2, median) * 100)), "points", sep = " ")
fd_changetext_short <- paste(fd_changelab_short, abs(round(apply(m1.fd, 2, median) * 100)), sep = "")

m1.fd_annotate <- data.frame(x = apply(m1.fd, 2, median), 
                             y = rev(unique(fd_long$Varpos)), 
                             out0 = paste(round(m1.fd_out0 * 100, digits = 1), "%", sep = ""), 
                             outROPE = paste(round(m1.fd_outROPE * 100, digits = 1), "%", sep = ""),
                             change = fd_changetext_short)
fd.p <- fd.p + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.1, size = 3) + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)

ggsave(fd.p, file = "Output_Tables-and-Figures/claims_m1_fd.pdf", width = 6, height = 8)

# Remove transparency
fd.p2 <- fd.p + annotate(geom = "rect", ymin = 0 - ROPE, ymax = 0 + ROPE, xmin = 0, xmax = Inf, fill = "gray90", color = NA) + geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = outROPE), color = "white", nudge_x = 0.1, size = 3)  + 
  geom_text(data = m1.fd_annotate, aes(y = x, x = y, label = change), color = "black", nudge_x = -0.1, size = 4)
ggsave(fd.p2, file = "Output_Tables-and-Figures/claims_m1_fd.eps", width = 6, height = 8)
